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Three-dimensional unsteady viscous effects can significantly influence 
the performance of fixed wing and rotary wing aircraft. These effects are 
important in particular for the flow about a helicopter rotor in forward 
flight and the flow about a three-dimensional (swept and tapered) 
supercritical wing at transonic speeds. In both cases, the viscous boundary 
layer is predominantly thin and exhibits regions of reverse flow in the 
streamwise and/or spanwise direction. Clearly a computational procedure that 
can calculate such flow fields would be of great value in the design process 
as well as in understanding the flow phenomena. 

A new computational procedure specifically designed to handle such 
problems has been developed which bridges the gap between the 
inviscid/boundary layer and Navier-Stokes approaches in that it is of 
sufficient generality to compute regions of reverse flow yet is efficient and 
user-oriented. Hence, the major objectives of the current effort were to 
adapt this technique to treat three-dimensional unsteady turbulent flows and 
validate the procedure by comparing the predictions to experimenal data. 

These goals have been achieved by considering the experimental data of 
Karlsson, i.e. two-dimensional unsteady oscillating turbulent flow over a 
flat plate. Two-dimensional calculations were performed and the results 
agreed both qualitatively and quantitatively with the data. Thereafter, the 
analagous three-dimensional case was considered which was obtained by a 
coordinate rotation to yield the flow over a flate plate skewed at 45° to the 
freestream direction. The results of this computation also agreed well with 
both the two-dimensional results and Karlsson f s data, hence validating the 
computational procedure in three dimensions. In addition, new inflow 
boundary conditions were developed and an explanation was proposed to resolve 
the discrepancy concerning other previously reported predictions of the skin 
friction phase lead angle as a function of reduced frequency. 

The computational procedure which has been validated in the Phase I 
effort can have a significant impact in the design of both rotary and fixed 
wings by predicting the onset of separation and the flow properties through 
regions of reverse flow. In addition, with the inclusion of a displacraent 
thickness interaction option, the procedure can be made more general and can 
aid in transonic flutter and buffet aerodynamic analyses where interaction is 
important. Furthermore, the calculation procedure can be extended to handle 
oscillating control surfaces in a three-dimensional flow. 
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ANALYSIS 


Introduction 

In recent years there has been increased attention given to three- 
dimensional unsteady aerodynamics. Such flows manifest themselves over fixed 
wing and rotary wing aircraft. Fixed wing aircraft are designed to be 
nominally steady with the unsteady effects being introduced through either 
control surface motions or induced external oscillations. These phenomena 
are observed in flow about wings throughout the Mach number regime; subsonic, 
transonic, and supersonic. The unsteady effects will influence loss levels 
as well as lift and moment coefficients. As the speed increases beyond 
subsonic into the transonic regime, additional effects such as flutter, 
buffeting and aileron buzz may arise. Since the new generation of energy 
efficient transports will be operating in the transonic range, the prediction 
of the onset of these phenomena under transonic conditions will play a 
significant role in the design process. Furthermore, the accurate prediction 
of the aerodynamic loading during the unsteady flow phase is crucial in 
prescribing the correct input for aeroelastic analyses. As a result of the 
three-dimensional nature of the transonic supercritical wings i.e. taper and 
sweep, the flow structure will be three-dimensional, while the oscillating 
recompression shock will cause regions of reverse flow to exist in the 
streamwise and/or spanwise directions. 

In contrast to fixed winged aircraft, the flow about a rotary wing is 
designed to operate in an unsteady environment. The flow about a helicopter 
rotor in forward flight is periodic and exhibits substantial unsteady 
transonic effects. As a consequence of the vortex-wake interaction of the 
advancing blades, broad band frequencies are excited having large amplitude 
oscillations relative to the blade. This flow is also highly three- 
dimensional with regions of reverse flow existing in the streamwise and 
spanwise directions. 

Concurrent with these developments has been an increased effort to 
better understand these phenomena by conducting three-dimensional 
unsteady wing experimental programs (cf. Ref. 1, 2) and applying inviscid 
computational procedures to these problems (e.g. Ref. 3, 4). The results of 
the numerical calculations in conjunction with the experimental data 
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indicate that the observed phenomena are strongly influenced by viscous 
effects near the body surface which are not accounted for by the inviscid 
predictions. These viscous effects are concentrated within a region that is 
predominately thin except for localized regions of reverse flow in the 
streamwise and/or spanwise directions. Hence, there is clearly a need to 
compute these viscous effects in an efficient and economical manner. 

There are several possible approaches available for computing 
three-dimensional viscous flows, ranging from empirical models to 
sophisticated treatments based on the solution of the three-dimensional 
time-dependent Navier-Stokes equations. Due to the complex structure of the 
flow the empirical approach is too restrictive. At the other end of the 
spectrum is the three-dimensional time-dependent Navier-Stokes analysis. 

Even though such procedures have been developed at SRA and have been applied 
successfully to a variety of problems (Refs. 5, 6, 7), such a technique is 
not required for many of the viscous layer type problems ocurring on wings in 
which the static pressure is sensibly constant across the viscous layer. 
Therefore, an approach is sought which allows the static pressure to be 
imposed at the boundary layer edge, but which can be used in 
three-dimensional flows having streamwise and/or crossflow separation. 

When the flow is steady and contains little or no separation, standard 
boundary layer prediction schemes have reached a high level of sophistication 
and predictive accuracy, even in three space dimensions. In unsteady flows, 
such as are commonly encountered in rotary winged aircraft and fixed winged 
aircraft at transonic speeds, some progress has been made in two space 
dimensions but little to date has appeared on unsteady three-dimensional 
boundary layers, especially for cases where negative cross flows are 
encountered, i.e., the spanwise component of velocity changes sign. To be of 
practical value, time-dependent three-dimensional boundary layer prediction 
schemes require high computational efficiency and transient accuracy coupled 
to the ability to treat arbitrary cross flow profiles. Conventional boundary 
layer integration schemes have developed by forward marching the streamwise 
velocity in the streamwise direction and simultaneously marching along the 
positive spanwise direction. However, in cases where there is "reverse cross 
flow", one must resort to specialized differencing to march a solution 
through such a region. 
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Recently, conventional boundary layer developers have applied an 
implicit spanwise construction to remove the restriction of only positive 
cross flows. Lin and Rubin (Ref. 8) in their predictor-corrector boundary 
region solutions for flow over a yawed cone at moderate incidence have also 
shown that retaining diffusion in the spanwise direction not only eliminates 
the problems associated with negative cross flows, but improves upon the 
solutions obtained by standard three-dimensional boundary layer techniques. 
Furthermore, with these spanwise implicit approach boundary conditions 
applied at the tip the flow inboard would be influenced, if required by the 
physics of the flow. 

As a consequence of these observations, a spanwise implicit formulation 
(retaining spanwise diffusion) is essential for both rotary and fixed wing 
applications, and it can be obtained for a very modest increase in code 
computational labor. Based on the experience in Refs. 9 and 10, the spanwise 
implicit sweep only results in a moderate increase in computational effort 
relative to the explicit spanwise marching approach. The extension of the 
conventional three-dimensional boundary layer equations to allow spanwise 
diffusion is easily accomplished, and in view of the improved physical 
representation, it has been implemented in this effort. 

Since the solution is being time marched, the opportunity to take a 
streamwise implicit sweep at roughly the same cost as the explicit sweep 
(forward march) does arise. If an implicit streamwise structure is adopted, 
then full time linearization can be utilized. That is the linearization of 
the nonlinear terms is performed about the known time level, rather than a 
known spatial marching level. As is pointed out in Ref. 10, it is easier to 
obtain a consistent high order accurate spatial-temporal linearization by 
marching in time than in space. Further, by structuring implicitly in the 
space marching direction, regions of axial reverse flow would not present a 
computational stability problem. As a result of these combined benefits an 
implicit structure in all three spatial directions has been employed in this 
effort. 

Since a full 3— D spatial integration is carried out at each time step of 
a transient calculation, spatial accuracy plays a very important role in the 
overall efficiency of the numerical method. The major goal of any spatial 
differencing scheme is to achieve a desired level of accuracy for the mimimum 
number of grid points. In addition, the scheme should be able to suppress 
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oscillations associated with coarse grid calculations, i.e. cell Reynolds 
number stability condition, without overly smearing the solution; a 
consequence of adding excessive artificial dissipation. Such a method (QR 
operator), to be discussed in subsequent sections, has been developed and is 
described in detail in Refs. 11 and 12. The QR operator scheme permits among 
its many options fourth order accurate spatial approximations. 

In this report the implementation of a computer code for the efficient 
solution of three-dimensional time-dependent viscous flows on fixed and 
rotary wing aircraft is described. The Linearized Block Implicit (LBI) 
technique of Briley and McDonald (Refs. 13 and 14) is employed in conjunction 
with QR operator technique to solve the present approximate form of the 
turbulent Navier-Stokes equations which are derived for nonorthogonal 
coordinates in generalized tensor form. The rationale for the choice of this 
approach is discussed in detail in Refs. 14 and 15. 

The basic assumption made in the derivation of these equations is that 
the pressure does not vary normal to the shear layer, and is obtained from an 
inviscid analysis. Inherent in this assumption is that the shear layer is 
thin. Since the boundary layer remains thin over most of the wing, except 
perhaps near the tip, it is expected that this assumption should not limit 
the applicability of the present approach. 

It is also assumed that for the energy equation the stagnation 
temperature, T 0 is constant. This assumption is a good approximation for 
the flow fields considered as discussed in Ref. 7, and is included here only 
for purposes of computer run economy. In the analysis that follows, the full 
energy equation could equally well have been used with consequent increase 
in computer run time. Employing the equation of state which relates the 
pressure p to the velocity components u and w by an algebraic equation, the 
problem with T Q assumed constant can be reduced to one involving only the 
three velocity components, u, w and v and three equations, the streamwise and 
spanwise momentum equations and the continuity equation. Hence, a 
block-three system is considered. If the full energy equation were to be 
considered, a block-four system would result due to the inclusion of the 
temperature as an additional unknown and thus would result in an increase in 
computer run time. 
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For turbulent flows, a two-layer mixing length model is employed and its 
formulation in generalized tensor notation is given. A novel method is 
employed for solving the continuity equation in conjunction with the momentum 
equations. In the following, a description of the computational procedure is 
given including, coordinate systems, governing equations, turbulence model, 
and numerical technique, i.e. QR operator and Linearized Block Implicit 
schemes. Thereafter, the computer code is described. Following this a 
detailed discussion is presented of the computations conducted, and the 
results obtained in meeting the Phase I objectives. 

Coordinate System 

Since the goal of this effort is to solve for the flow over wings and 
rotors an understanding of the type of geometries to be considered is 
essential to guide the choice of the coordinate system and the structure of 
the computer code. The coordinate system is not only dependent upon the 
geometry of the wing, but also upon the approximations that are made to the 
governing Navier-Stokes equations. As in boundary layer theory, the present 
approach assumes the pressure is constant normal to the shear layer. 

Inherent in the assumptions is that the shear layer is thin. As pointed out 
by Howarth (Ref. 16) the boundary layer assumptions lead to the conditions 
that one coordinate direction must be normal to the body surface while the 
other coordinate directions must lie on the body surface. These conditions 
uncouple the metric data on the surface from that in the normal direction. 
Hence, the metric data for the surface coordinates are functions of the 
surface coordinates alone, while the metric data for the normal coordinate 
direction are functions of that coordinate alone. The choice of the surface 
coordinates is rather arbitrary and is based on considerations such as the 
ease of construction of the grid distribution on the wing surface. 

For example, a rectangular planform wing could be adequately described with a 
Cartesian coordinate system. However, for more general planforms such as a 
swept and tapered wing a nonorthogonal grid which conforms with the 
boundaries is preferred since it represents the airfoil more accurately 
(cf. Fig. 1). 

Another consideration is the selection of a coordinate grid 
distribution; the major objective being the resolution of large solution 
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gradients. The approach taken here is to construct coordinate trans- 
formations that contain distributions for physical mesh points. In this 
context, the uniform mesh of computational space is simply mapped into a 
suitable distributed mesh in physical space. When the transformation 
contains the mesh point distribution, there is no need to construct the 
apparatus for the discrete approximation of derivatives on a nonuniform 
mesh. This results in a savings in both computer logic and storage. 

Hence, in this work a coordinate system is chosen that conforms with 
the boundaries of the physical domain, i.e., the wing surface which in 
general will be nonorthogonal with provisions made for analytical grid 
transformations (Ref. 17) in each coordinate direction. 

In view of the type of geometries to be considered and the assumptions 
made to obtain the equations, a specialized nonorthogonal coordinate system 
is advocated where the metric tensor which has four independent components is 
given by 

r Qii 9i2 o i 


J ij 


9.2 9 22 


o 0 9 33 


The subscripts 1 and 2 refer to the directions on the surface of the body 
while subscript 3 refers to the direction normal to the body. Since the 
metric data in the coordinate directions on the airfoil surface are not 
functions of the normal direction, the metric data in a parallel surface 
above the body are evaluated on the body surface (Ref. 16). Furthermore, due 
to the use of nonorthogonal coordinates it becomes advantageous to derive the 
equations in general nonorthogonal coordinates employing generalized 
tensors. Details of the generalized tensor notation can be found in 
Refs. 15, 18, and 19. 

An important feature of the analysis to follow is that the governing 
equations which are derived, under the prescribed assumptions, are invariant 
for any coordinate system or any grid transformation (although, of course, 
the physical approximations are coordinate dependent). The grid 
transformations are absorbed into the geometrical coefficients, leaving the 
equations unaltered in form. This feature allows for the construction of the 
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geometric data to be contained in a single subroutine with the definition of 
the metric data and their derivatives as input. 

In Ref. 21 the geometric properties of a surface in three-dimensions 
are discussed and where appropriate, the generalized tensor equivalents are 
given. In addition, a symmetric (uncambered) NACA four-digit airfoil is 
considered and the pertinent geometric coefficients are presented. 

Governing Equations 

In the following, the governing equations are nondimensionalized 

as follows, x* with respect to the characteristic length L, the velocity 

with respect to Uoo, density, pressure and temperature with respect to po>, 

2 2 

p OT U M and Uoo /cp respectively and time with respect to L/Uoo. The 
viscosity is nondimensionalized with respect to Uoo. 


Continuity Equation 


dp 


at 


+ 


j 
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where J is the Jacobian, p the density, and u^ is the kth contravar iant 

velocity component 
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where ' ,k' denotes a partial derivative, ' denotes a covariant 
derivative and g^ is a component of the metric tensor and A is the 
velocity divergence. 

In Ref. 21 it was pointed out that the QR Operator scheme requires that 
the governing equations be in quasi-linear form and that the spatial operator 
in a given direction operate on only one variable. For the momentum equation 
this requirement prevents the implicit treatment of certain diffusion terms 
that arise due to the curvature effects. In the usual boundary layer 
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approximations these explicitly treated terms would not appear in the 
equation since they are of order 0 (Re“^/^) or smaller, and should, 
therefore, be of little consequence. 

Since mixed partial derivatives are commonly treated explicitly in 
orthogonal coordinate systems, we do likewise in generalized nonorthogonal 
coordinates and extend this concept to include mixed second covariant 
derivatives. All other second covariant derivatives are retained as 
implicit. Since the pressure is specified and impressed upon the viscous 
layer, its specification replaces the normal momentum equation. Thus, the 
streamwise and spanwise momentum equations are the only two retained. A more 
detailed discussion of the derivation of these equations is given in Ref. 21. 


Energy Equation 


For the energy equation constant stagnation temperature is assumed. 
Neglecting the square of the normal velocity with respect to the squares of 
the other velocity components 


T = T + 


t< U P 


2 + w *) + 


'12 


h,h 2 


u_w_ 

P P 


(3) 


where Up and Wp and the physical velocity components. These assumptions 
are employed here only for simplification purposes. If warranted, they can 
be removed and the full energy equation can be considered. 


Equation of State 

The equation of state assumes a perfect gas and is given by 

r - 1 


p = 


7>T 


(4) 


Linearizations 

The following analyses assume a set of linear partial differential 
equations. However, the convective part of the momentum equation and the 
continuity equation are nonlinear, containing terras that involve the product 
of density and velocity components. In order to overcome this difficulty, 
the procedure described in Ref. 10 is employed to linearize the 
aforementioned terms by Taylor series expansion about the known time level 
solut ion. 
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It is important to note that in the governing equations the 
contravar iant velocity components are used. However, as noted in Ref. 22, it 
is advantageous to solve for the physical velocity components. Therefore, 
when the governing equations are subsequently cast into a form amenable to 
the application of the LBI scheme, they are transformed so that the physical 
velocity components appear. 

Turbulence Model 

In turbulent flow cases, the three-dimensional ensemble-averaged 
turbulent flow equations are considered. The approach taken here assumes an 
isotropic turbulent viscosity, p^, relating the Reynolds* stress tensor to 
mean flow gradients. Using Favre averaging (Ref. 23) the governing equations 
then are identical to the laminar equations with velocity and density being 
taken as mean variables and viscosity being taken as the sura of the molecular 
viscosity, p, and the turbulent viscosity, p-p. 

At this point additional closure assumptions for the Reynolds stresses 
are required, i.e., the evaluation of Pj. There are a variety of 
approaches available, from the simpler mixing length models to the more 
complicated one and two-equation models. Since the intention here is to 
verify the code’s performance in wall bounded cases, the mixing length model 
which has worked well in the past for similar flow environments (Ref. 24) was 
chosen* The extension to more complex models could be undertaken at a later 
time if warranted. At that time, the LBI procedure that is used for the 
solution of the momentum equation could be applied to the kinetic energy of 
turbulence and the . turbulent dissipation equations. 

Employing the Prandtl mixing length concept, the turbulent viscosity is 
given as 


m ■- pi 2 y$ w 

where l is the mixing length and <t> is the dissipation function, which in 
generalized tensor notation is given by 
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As in the Cartesian tensor formulation, $ does not automatically reduce to 
the dominant term for standard boundary layers, i.e., 3u/3y in two dimensions 
and ((3u/3y) 2 + (3w/3y J 2 ) 1 / 2 in three dimensions. Hence, provisions are 
made in the computer code that on option retains only the dominant components 
of the strain. 

The mixing length formulation is based on McDonald's model (Ref. 25), 
and is given by 

l -- I^q nh % 

*00 ( 8 ) 

where £« is the outer layer length scale and 

= I - exp(-y + /A + ) (9) 

where y + takes on its usual meaning. The constants appearing in Eqs . 8 and 
9, k,A and A + are .4, .09 and 26.0 respectively and 6 is the local boundary 
layer thickness defined as .995 U e . Note that in the limit as y+0 Eq. 8 
reduces to 

/,= ky3 

while for y large Eq. 8 reduces to 



the standard two layer values. 

Spatial Difference Approximations 

QR Operator Notation 

In this section, implicit tridiagonal finite difference approximations 
to the first and secdnd derivatives and to the spatial differential operator 
are considered. The QR operator procedure for generating a variety of 
spatial discretizations is also introduced. As special cases, standard 
second-order finite differences, first-order upwind differences, fourth-order 
operator compact implicit (OCI), fourth-order generalized OCI and exponential 
type methods are obtained. Since all these schemes are of the same form 
(cf. below), a single subroutine which defines the difference weights is all 
that is required to identify the method, while leaving the basic structure of 
the program unaltered. The rationale for the use of the QR approach in the 
present problem is discussed in detail in Ref. 14. 
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The QR formulation allows for ADI methods and permits the treatment of 
systems of coupled equations, i.e., LBI methods. Although variable mesh 
schemes can be employed within the QR framework, it is believed preferable to 
use analytic transformations to obtain a uniform computational mesh, hence 
attention is restricted to uniform mesh formulations. 

The general concepts and notation will be introduced for two-point 
boundary value problems and then the methodology will be extended to more 
general linear and nonlinear parabolic partial differential equations in one 
dimension. The application of QR operator method to multidimensional 
problems is discussed in the section pertaining to the LBI scheme. 

Consider the two-point boundary value problem 

L(u) * a(x)u xx + b(x)u x + c(x)u * f(x) do) 

with boundary values u(0) and u(l) prescribed. Derivative boundary 
conditions, although not discussed here, can easily be incorporated into the 
framework of the QR operator notation. Let the domain be discretized so that 
xj = (j-l)h, j= 1, 2, . . . , J + 1, and Uj~u(xj), F j ~ 
u x (xj), Sj — ^u xx (xj) and h = 1/J is the mesh width. The numbering 
convention was chosen here to be compatible with FORTRAN coding. 

Without loss in generality for a(x) * 0, Eq. (10) can be divided by 
a(x) so that we may treat instead the following equation 

L(u) » u xx + b(x)u x + c(x)u»f(x) (11) 

where 

b(x) - b(x)/a(x), c( x ) - c(x) /o(x) and f(x) - f(x)/a(x) 


Substituting the finite difference approximations to the first and 
second derivatives 


U i+ .- u i . ? 

8iT“J ' 2 h =F i- U « (X i ) + 0lh > 


u. 


D + D_ 


U: 


Uj-r^Uj+Uj+i 


h 2 h 2 

into Eq. (11) and rearranging, we obtain 


• s J -u xx (x j ) + o(h 2 ) 


( 12 ) 

(13) 
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L(u) - 



or 



fj (14) 


[ 


I - 




(15) 


where Rcj = hbj is the cell Reynolds number. 

Equation (15) can be generalized by introducing operator format, i.e., 

rj Uj., + r J Uj + rj U j4l - h 2 (q]fj_, + qjfj + qjf J+l ) . (16) 


where the superscripts (-) minus, (c) center, and (+) plus indicate the 
difference weight that multiplies the variable evaluated at the (j-1), (j) 
and (J+l) grid points respectively, and where the rj's and qj's for grid 
point j are functions of h, bj_j, b j , bj+i, cj_i, Cj and cj + j. 

Comparing Eqs . (15) and (16) we can identify the rj's an d *1 j * s » viz.. 


rj - I - RCj /2 


qj ■ O 


r j * ^ 2c J ~ ^ 


fj - I + Rcj /2 


"j ■ 1 

qj ■ 0 

We now define the tridiagonal difference operators Q and R 

r[uj] ■ r j Uj_, + 'j Uj + r ; Uj„ 

o[fj]"qjfj-i +£ ij f j + Rj f j«i 


(17) 


(18) 


Noting that L(u) = f and substituting Eq. (18) into Eq. (16), we obtain 

R [ U j] -h 2 o[ f j]-= hS °[ l_(u)j] (19) 
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Alternatively by employing the inverse operator an expression for 

L(u)j can be obtained 


L(u)j * 



0' ! RUj 


( 20 ) 


For standard central finite differences Q = Q“* = I, the identity 
matrix, so the spatial operator can be given explicitly in terms of Uj_i, 

Uj and Uj+i. However, in general, for higher order methods whereas Q is 
tridiagonal Q -1 is a full matrix, and the spatial operator cannot be given 
explicitly in terms of the variables at adjacent grid points. Hence, Eq. 
(20) provides a method for expressing the spatial operator for a wider class 
of difference approximations. The formalism in Eq. (20) is also applicable 
for first and second derivatives appearing alone (cf. Ref. 11). In Refs. 12 
and 14 a technique due to Berger, et al is described for constructing fourth 
order tridiagonal methods which possess a raonotonicity property as the cell 
Reynolds number is increased, R c + 00 . This type of scheme is an option in 
the computer code. 


Application to Coupled Nonlinear Parabolic Equations 

Before considering the LBI technique, we discuss some of the 
limitations placed on the QR operator scheme in solving a system of nonlinear 
parabolic equations. 

Given a system of m nonlinear parabolic equations in ra unknowns, 


m 


i 


°r /J 


. n+i n. 

< u ij - u ij> 

At 


n+/3 # 


m’ 


'j > *2 9 ^ ^ 


0 


n+B 

where is a 
directly over 
operated upon 


j - 1,2, J+l 

quasilinear spatial operator, the QR formalism carries 
provided that for any equation only one independent variable is 
by the differential operator. For example, 


I 

a(u,w,v) u t 


+ b(u,v,w)u x + c(u,v,w) 
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is allowed since x derivatives of u only appear, while 

— 7 “T u l * u vv + b(u,V,w)u y + c(u,V,w) + dCUjVjWjWj- 

atUjWjV) 1 xx x * 

is not allowed since x derivatives of both u and w appear. The approximate 
form of the unsteady Navier-Stokes equations used here, when written in 
quasi-linear form, falls within the class of allowable differential 
operators. Thus, for the problem being addressed in the present study, the 
OCI schemes are applicable. Note that within the splitting approach, 
nonallowable terms in the OCI scheme such as dw x above, may be split off 
and treated by a special implicit sweep. Provided care is taken and for 
instance the Douglas-Gunn (Ref. 26) formalism is adhered to, no particular 
problem arises other than the cost of an additional implicit sweep which is 
incurred. 

Thus, multidimensional problems and/or more general equation forms 
can usually be accommodated by a splitting procedure, which reduces the 
differential operator to a sequence of one-dimensional problems which 
have the appropriate allowable form. However, as with standard finite 
differences, to avoid the cost of additional implicit sweeps, special 
procedures must be applied to cross derivative terms, e.g., extrapolation or 
explicit treatment . 


Linearized Block Implicit Scheme 
Consider a system of nonlinear partial differential equations 

A$ t = (21) 

where $. is a vector of unknowns and ? is a source term vector which is a 
function of x 1 , x 2 , x 3 and t. Extension to source terms which are functions 
of $ are discussed in Ref. (10). is a three-dimensional nonlinear 
differential operator and the matrix A appearing in the momentum equations is 
equal to pi where p is the density and I the unity matrix. 

Equation (21) may be centered about the n+B time level, i.e. t n+ ^ - 
(n+B)At = nAt+BAt = t n +BAt, and written 


A n+0[3!* , -5 n ] / At 


n+/3-rn+/3 -n+£ 


( 22 ) 
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where 0 _< M 1 is a parameter allowing one to center the time step, i.e., 
0=0 corresponds to a forward difference, 0 = 1/2 to Crank-Nicolson and 0 = 

1 to a backard difference. 

After linearizing Eq. (22) by Taylor series expansion in time about the 
time level by the procedure described in Ref. 10 to give a second-order 
linearization, we obtain 

A n [5" t -e' n ]/At = ^ n [$" t A.$"]-2, n * n + ( 23 ) 

where X is the linearized differential operator obtained from lb. 

The difference between the nonlinear operator and the linear 
operator 1 is defined as - l n . At the intermediate level 

n + B, $n+B £ s represented as 

= ^5 n+ ' + (1-/3)$" (24) 

Using these relationships and dropping the vector superbar for convenience a 
two-level hybrid implicit-explicit scheme is obtained 

A n (4> n+ '-<J) n )/At = ^"«>" + '-4>") +^"4> n + m" <t>" + (25) 

The vector \|> n+ P represents all of the terras in the system of 
equations which are treated explicitly. More about this will be said later, 
but for the moment note that ma y be approximated to the requisite 

order of accuracy by some multilevel linear explicit relationship, or 
approximated by with a consequent order reduction in temporal accuracy. 

The operator X is now expressed as a sum of convenient, easily 
invertible suboperators 1 = 2 \ + 1 2 + . . . . J m . In the usual ADI 
framework these suboperators are associated with a specific coordinate 
direction. Further, it is supposed that these suboperators can be expressed 
in the QR notation introduced earlier. Writing iJ> n+ S anc j M n $ n as a 
single source term S n+ &, Eq. (25) is written as 

a"[ /a« +^ 3 n ][<j >n+ -<j> n ] + [j^j n ^r 3 ]i t> n + s nt/3 (26 ) 
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To solve this system efficiently it is split into a sequence of easily 
invertible operations following a generalization of the procedure of Douglas 
and Gunn (Ref. 26) in its natural extension to systems of partial 
differential equations. The Douglas-Gunn splitting of Eq. (26) is written as 
the following three-step procedure 

A "[ $♦-$"] /ai + [ J ?, n + i 2 n +4 r, ]4>"+ 

A n [4>*'-$ n ]/it = /<[<*>* -<i + /3J "[•*•"- <*>"] + ik + s "^ 
A " [$”-4."]/" = /St/fl®* -<*>"] + /9^[4”'-4 n ] 

+ [/" "+^ 3 n ]<t>" + S n *& (27) 

which can be transformed to the alternative form 

f n 

[ A 


-AtW,"] [$*-<*>"] ■ Al[/%^ 2 "^ 3 n ]<J. n + AtS n ^ 
[a" - At/3. 1 ”] [4.*“-4. n ] = a"[ 3>*- <3> n ] 

[a"- At /3^ 3 n ] = A n [<t> ,M, -4 n ] 


If the intermediate levels are eliminated, the scheme can be written in the 
so-called factored form 

(A"-PAl/ 1 "XA n )''(A ,, -/3At/ 2 n )(A")’'(A n -)3At/ 3 n )(<J> nM -$") = (29) 

At( J? + jf 2 " +/ 3 ") 4>" + At S n * /3 

The ADI formulation given in Eq. (28) is directly applicable for J! ^ 
operators represented in Q“*R operator format. Consideration of 
intermediate boundary conditions and the removal of the inverse operator 
Q”1 is given in Ref. 21. 

It is worth noting that the operator or J! can be split into any 
number of components which need not be associated with a particular 
coordinate direction. As pointed out by Douglas and Gunn (Ref. 26), the 
criterion for identifying sub-operators is that the associated matrices be 
’’easily solved" (i.e., narrow-banded). Thus, mixed derivatives and 
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complicating terras which might inhibit the use of OCI can be treated 
implicitly within such a framework, although this would increase the number 
of intermediate steps and thereby complicate the solution procedure. 

An inspection of Eq . (28) reveals that only the linearized operators 

Jt. i 2 and 2 appear. Indeed, the computer code employs this feature by 
1 2 3 

evaluating these three operators before the first sweep, storing them and 
accessing them as needed in the subsequent three sweeps. In addition, the 
terms arising from the nonlinear terms are immediately absorbed into S n+ ^ 
as they appear, allowing for an efficient evaluation of the terms in the 
differential equations. 


The spatial operators appearing in the differential equations 

•^1* 2 and 3 must be identif ied at least formally in order to isolate 

the coefficients that are to be used in the construction of the Q and R 
operators. These operators can be represented in standard form at each grid 
point, i.e.. 


J> n = a n $ + a n <I> 


i i 


i,ii 


12 




+ a n + o n + a n <J> 

13 I 14 2 U |5 - 


(30) 


In Eq. (30) the first subscript of 0 indicates the velocity component 
(associated with the corresponding direction) and " , " indicates a 

derivative. The subscripts of the a”^ refer to the direction (i) and the 
term in the equation (j) respectively. Note that the equation is in 
quasi-1 inear form, since the coefficients of the derivative operators need to 
be identified, for use with the QR operator technique employed here. 

Alternate schemes have been proposed by Leventhal (Ref. 27) for equations in 
conservation form but are not considered here. In the following section, a 
description will be given of how this entire operator is discretized by 
employing the QR operator format, and how the discretization is incorporated 
into the LBI framework in order to solve the system of equations (28). 

The continuity equation is considered first. Since it is a first-order 
partial differential equation it does not have the standard form of Eq. 

(29). Furthermore, in the linearization process p has been eliminated in 
favor of the u 1 velocity components so that the continuity equation has 
become an equation for the three velocity components, and not density. 
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An inspection of the system of equations under consideration reveals 
that substantial savings can be realized if the equations are partioned 
appropriately. Since the normal to the body ( M 3 M direction) momentum 
equation is not solved, the normal velocity appears only in conjunction with 
terms associated with the normal # '3" direction in the other two momentum 
equations. Hence, in the first two sweeps where directions ”1" and “*2 H are 
implicit one is required to solve only for the two corresponding velocity 
components in the streamwise and spanwise momentum equation without the need 
of considering the continuityequation. However, on the third sweep where all 
3 velocity components appear, one must solve all 3 equations. This strategy 
reduces the solution procedure to the inversion of two 2x2 block matrices 
and one 3x3 block matrix rather than three 3x3 block matrices which leads 
to a substantial reduction in computation time. If the full Navier-Stokes 
equations were considered (including a normal momentum equation) the 
aforementioned partioning could not be applied since the normal velocity 
would appear in all three sweeps. 

The question that arises is how to appropriately split the continuity 
equation, since it need only be solved on the third sweep. Here again the 
Douglas-Gunn formulation leads to the appropriate choice. The continuity 
equations written in conservation form is, 



After linearizing and eliminating p, the increment form is obtained 

A n a n+i , n n Au, n+| , At/3 d r n A n A n+i , n^n A n+i, n A n+fl 

A Au + B AW + — — ‘ ^ 3 - [v A Au + v B Aw + p Av J 

At d r i I In . At)3 d T / n j. n a°\ a n+i, . n_n. A n+H 

= L J / >u J + — ~dJT L { / > + u A )Au +(u B )Aw J (32) 

d r, n , n„n x A n+i , , n 4 n %A _n+r| 

+ — - — [(/> + w B )Aw +{w A )Au J 

where all the velocity components are the contravariant components u = u 1 , 

2 3 

w = u and v = v . J is the Jacobian and 
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A " = K/ + 9,2 W "] 
b" = |;r[g22w n+ g,2u"] 


By employing the Douglas-Gunn procedure, Eq. (32) is represented as a 

third sweep equation, and a consistent approximation is obtained to the 

continuity equation, i.e., the x 1 derivative term is evaluated at the * level 
2 

and the x derivative term is evaluated at the ** level. The values of the 

intermediate derivative terms are obtained after the solution of the first 

two sweeps of the two momentum equations. Note that these terms do not 
contain the normal velocity* The equation can thus be written in symbolic 
form 

n A n+i n n+i At/3 d r r n n A n+i n n . n+i , n . n+i 1 1 

A Au + B Aw + — — [ J(A v Au + v B Aw + p Av | J 

X (33) 

= s "-^4r [*{}}* -&[*{}]” 


• • o 

Since the only terra involving v is in the x derivative terra, one can 
directly integrate the equation with respect to x , i.e. 



The next section describes how this is done very easily via the QR operator 
scheme. The concept of integrating directly the continuity equation is not 
new. Blottner (Ref. 28) attributes a similar coupled procedure to Davis for 
the two-dimensional steady boundary layer equations in which the trapezoidal 
rule is used to integrate the continuity equation. Weinberg (Refs. 29 and 
30) also used a fourth-order Simpson integration scheme to solve the compres- 
sible boundary layer equations. Such procedures are stable and offer a 
viable alternative to approximating the derivatives by finite differences. 
Note that conceptually the continuity equation in integrated form is treated 
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on each sweep of the Douglas~Gunn splitting, although in actuality the 
continuity equation is solved only on the third sweep* The stability and 
consistency of the original splitting is still retained since the split 
integration operators can be incorporated into the J and difference 
operators (cf. Eqs. 27 and 28). 

Implementation of the LB I Scheme Employing 
the QR Operator Technique 

Consider the third sweep of Eq. (27) in which both momentum equations 
and the continuity equation are solved. The momentum equations are in the 
form 

[a - /3At/ 3 JA4> = A A<J> (34) 

where A$*** is the column vector of unknowns, u, v, w* Here it has been 
implicitly assumed that the equations have been appropriately normalized and 
that the contravariant velocity components have been suitably transformed 
into their physical components. Employing physical components, (cf. Ref. 22) 
leads to a better behaved solution since these components are not unduly 
influenced by geometrical variations. 

For the streamwise momentum equation one obtains 

A <£***= Au, 33 + u 23 Au, 3 + a 33 Au + q 43 Aw + a 53 Av (35a) 

while for the spanwise momentum equation one obtains 

J" A<t>***= Aw, 33 + bjjAw.j + b 33 Aw + b„Av + b 5 ,Au (35b) 

where superscript *** has been omitted from Au, Av and Aw. Now in Eq. (35a), 
the first three terms on the right-hand side are approximated by the operator 
equivalent so that 

A" A®'”*. +0„Aw + 0 53 Av (36) 

3 
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Similar approximations are made for Eq. (35b). After substituting Eq. (36) 
into Eq. (34), and multiplying through by Q, one obtains for the streamwise 
momentum equation 

[Q,/» n - /?XR,]au - = Q,^> n Au** (37) 

2 

where A = At/Ax 3 . A similar expression is obtained for the spanwise 
momentum equation. 

The same type of procedure is also employed for the continuity 
equation. Since the continuity equation involves only first derivatives, 
they can be represented as 

d - QcRc (38) 

dx 3 Ax 3 

were Qc and Rc are the QR operators associated with the continuity equation. 
The operators Qc and Rc are constructed to approximate the weights associated 

with either a second-order trapezoidal rule or a fourth-order Simpson’s rule 
i.e.. 

Trapezoidal rule 



Simpson's rule 
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The discretized continuity equation thus becomes after multiplying thru by 
Qc 


[Q c JA n + £u)R c JA n v n ] Au + [o c JB n + )9tuR c JB n v n ]Aw + [ AV = Q C (RHS) (39) 

where RHS contains all the terms due to the linearization procedure and the 
terms evaluated at the * and ** levels and m = At/ Ax 3 * The resulting matrix 
derived from Eqs. (37) and (39) becomes a block 3 tridiagonal matrix 
(Ql , R i, Q2t r 2 » Q c and Rc are tridiagonal operators) with each sub block 
taking the form 


[ 0,p n -£Xr, ][ -/3AIO,a 43 ] 

[><**] 


Au 


Q,(Au**) 

[ -/?At Q 2 b 53 J £ Q 2 p j 

[-/3At0 2 b„] 


Aw 

= 

Q 2 (Aw**) 

[Q c JA n + /3cu R c JA n v n ] [Q c JB n +/3tuR c JB n v] [/3iuR c J/> n ] 

- V _ 


Av 


0 (RHS) 
c 


The matrix is inverted by standard LU decomposition. 

The Computer Code 

The type of numerical algorithm employed as well as its formulation has 
a marked impact on the structure of the computer code. One needs to consider 
both the number of CPU operations as well as the memory requirements. 

Usually, the number of operations can be reduced at the expense of increasing 
the amount of storage. For the type of problems under consideration both in 
two and three dimensions, with their large data bases, core requirements be- 
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come a stringent constraint and external disk storage must be employed. An 
additional consideration is the allocation and acquisition of data, in parti- 
cular how this impacts upon code vectorization. With these requirements in 
mind the code has been structured in such a way as to exploit available core 
and disk files by employing concurrent data transfer (buffer type statements) 
and structuring the code to operate on lines of data so that vectorization is 
possible. The calculations that are discussed in the next section were 
obtained on the CRAY 1 computer situated at the University of Minnesota 
computer center. The original version of the code was developed to run ef- 
ficiently on a CDC 7600 scaler machine. That version was transferred to the 
CRAY, but was not modified to exploit the vectorization capabilities of the 
machine since that was not within the scope of the present contract. Hence, 
the quoted run times, which are low, could still be reduced substantially if 
full vectorization is taken into account. This could be part of a Phase II 
effort. However, certain features have been incorporated into the code which 
leads to the mentioned efficient run times. The modifications described 
above are suitable for other vector processers such as the CYBER 205 as well 
as the CRAY. 

An investigation of the operation count of the LBI scheme in 
conjunction with the QR operator technique reveals that the most significant 
fraction of time is spent in computing the matrix coefficients, i.e. the 
linearization coefficients and difference weights, and exceeds the time 
required for the matrix inversions. Hence, it is worthwhile to optimize the 
calculation of these coefficients, and if possible store their values. This 
was accomplished by storing the operator coefficients and ./ as they 

were computed in the first sweep on the right-hand side of the differential 
equation. On the second and third sweeps X ^ and X ^ were accessed 
respectively and were not recomputed. It was for this reason that the 
formulation of the LBI scheme referred to the linearized operators J e 
instead of the ^>?'s on the right-hand side of the equation. 

The general structure of the computer code will now be described. 

After the input section and the initialization of data e.g. geometry, grid 
transformations, initial flowfield, etc. the actual construction of the 
difference operators is begun. The first derivatives of the velocity 
components and viscosity are obtained for the entire flow field and stored 
for ready access when needed for the computation of the appropriate terms in 
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the governing equations. Thereafter the terms that are to be treated 
explicitly are evaluated and absorbed into the function S n . 

The operators J! j, J, 2 and ^3 are then computed. These are 

used to evaluate the appropriate Q and R coefficients which are then stored 
for easy retrieval during each of the ADI sweeps. 

In the first sweep the matrix resulting from the application of the 
J! Y operators for the streamwise and spanwise momentum equation is 
solved as a 2 x 2 coupled system. The solution of this system, the * level 
quantities, are then used to construct the right-hand side of the second 
sweep equations and to evaluate the appropriate * level terra in the 
continuity equation. At this point the J! 2 operator is accessed and again 
a 2 x 2 system of equations for the streamwise and spanwise momentum equation 
is solved. The ** level quantities are then used to construct the right-hand 
side of the third sweep equations as well as the appropriate terms in the 
continuity equation. For the third sweep equations which consist of the two 
momentum equations and the continuity equation, the J! 3 operator is 
accessed from memory. The resulting 3 x 3 system of equations is solved for 
the three velocity components. 

After the primary variables are evaluated, the thermodynamic 
quantities, density, temperature and viscosity are computed. The procedure 
is then repeated at the following time steps. 

Discussion of Results 

In order to validate the numerical procedure for solving the 
three-dimensional unsteady turbulent approximate form of the Navier-Stokes 
equations it is necessaary to compare with available experimental data. 
Unfortunately, reliable three-dimensional data is rare and it is necessary to 
resort to alternate strategies in order to achieve this goal. Previously one 
of the present authors has successfully employed a procedure for verifying 
the numerical method for three-dimensional steady laminar flow (cf Ref. 14). 

A three-dimensional flow field was derived from a two-dimensional one by 
performing a coordinate rotation of the computational domain through a 
specified angle (for instance 45°) as shown in Fig. 2. As can be seen in 
Fig. 2 this is analogous to considering a uniform flow over a flat plate 
skewed at the rotation angle. In the orthogonal coordinate system consisting 
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of the leading edge of the plate and the normal to it (5,0 the flow is 
two-dimensional. However, when viewed in the x 1 , x 2 coordinate system the 
flow is a restricted type of three-dimensional flow involving in addition to 
the streamwise and normal velocities a spanwise component. Although this 
represents a restricted three-dimensional flow, it serves as a valid 
three-dimensional test case for the code. 

Such a procedure has several attractions over considering an actual 
three-dimensional flow. First, the analogous two-dimensional calculation (in 
the 5,5 coordinate system) can be performed giving additional information 
with which verifications can be made to assess proper code operation in the 
thee-d imens ional mode. Second, three-dimensional geometrical effects are 
removed allowing for a more direct comparison of the numerical results with 
the data. Third, turbulence modeling is simplified since the models have, in 
effect, no three-dimensional complications other than the contribution to the 
dissipation appearing in the mixing length formulation (cf Eq. 6). Finally, 
the two-dimensional calculation can facilitate the validation of the 
three-dimensional calculation. 

Having resolved the three-dimensional problem to be considered, the 
remaining task reduces to determining the appropriate two-dimensional 
experimental case with which comparisons are to be made. In recent years 
there has been great interest in the experiments of Karlsson (Ref. 31). 
Karlsson conducted a set of experiments for a zero pressure gradient flat 
plate flow whose external velocity consisted of sinusoidal fluctuations 
superimposed upon a mean velocity, i.e., 


U = U 0 ( I + A cos cut ) 

(40) 

where A is the amplitude of the oscillations and o> is its circular 
frequency. The magnitude of the amplitude ranged from 8% to 34%, while the 
oscillation frequency ranged from 0 to 48 Hz. Karlsson measured the mean 
velocity, the in phase and out of phase components of the first harmonic of 
the periodic fluctuations as well as the sum of the turbulent intensity and 
the contribution of the higher harmonics. A schematic of this test facility 
reproduced from a report by Carr (Ref. 32) is shown in Fig. 3 , and consists 
of a rectangular tunnel with the unsteady fluctuations being introduced by a 
rotating shutter assembly at the exit of the tunnel. 
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In Fig. 4 are shown recent computational results that have been ob- 
tained for the phase angle lead of the skin friction (cf. Refs. 32, 33, 34, 
35, and 36). The phase angle of the first harmonic of the skin friction is 
plotted against the reduced frequency, 0)x/u o . It is immediately apparent 
that there is a wide variation in the predicted results. For laminar flow 
there is a direct correlation between the phase angle and the reduced 
frequency. As can be seen in Fig. 5 (cf. Ref. 21), the phase angle ap- 
proaches an asymptote of 45° corresponding to Lighthill’s high frequency 
limit (Ref. 37) for large values of the reduced frequency. For turbulent 
flow, as shown in Fig. 4, there does not appear to be this direct correlation 
nor does the data approach an asymptotic value as observed for the laminar 
case. Such an asymptotic behavior would be plausible from physical 
consideration since the skin friction is dependent on the local flow 
properties near the wall, i.e. viscous sublayer. Furthermore, there does not 
appear to be any correspondence among the various predictions. An 
explanation for some of these disparaties is proposed and will be discussed 
below. For the moment it must be kept in mind that the experiments conducted 
by Karlsson are more than 25 years old and the methods used for obtaining and 
reducing the data were not as sophisticated nor as reliable as they are 
today. More disturbing from a computational point of view is that there is 
data at only one strearawise location and there is no upstream data with which 
to make comparisons. Another concern is the low Reynolds number, 
approximately 10 /ft at which the experiment was conducted. However, the 
the experimental data remains valuable and can be used to verify calculation 
procedures . 

In calculating this flow it was intended to be as faithful as possible 
to the actual flow conditions of the experiment. This could not be totally 
achieved due to the low Reynolds number of the experiments which would neces- 
sitate the consideration of unsteady transition, and the inclusion of low 
Reynolds number effects. Since this was not within the scope of the present 
effort, fully turbulent conditions were assumed. Hence, at the upstream 
boundary the velocity profile employed was a fully turbulent one. This is in 
contrast to some other numerical predictive methods in which it is assumed 
that the laminar flow developing from the leading edge of a flat plate 
undergoes instantaneous transition to turbulent flow at some predetermined 
strearawise location. 
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In view of the lack of upstream information and the low Reynolds 
number, the determination of a physical streamwise location for the upstream 
boundary required some attention. The experimental data of Wieghardt 
(Ref. 38) in conjunction with the curve fits of Clauser (Ref. 39) were the 
starting points for determining the upstream profile as well as the physical 
streamwise extents of the computational domain. In Fig. 6 Wieghardt 's data 
is plotted in Cf, Reg* coordinates along with Clauser 's predictions. As 
can readily be seen the experimental data compares well with Clauser 's 
•predictions for Reg* > 2000, but diverges from Clauser 's prediction at 
lower values of Reg* where low Reynolds number effects are prevalent. For 
the present calculations, the upstream boundary was placed sufficiently far 
upstream of the measuring station (nominally Reg* = 3600) where the 
numerical predictions are to be compared to the experimental data. This 
location was chosen to be approximately Reg*= 850, with the corresponding 
value of Cf = .00051. In Fig. 6 this point lies between the Clauser curve 
and Wieghardt 's data. Modifying these values by as much as 10% had little 
effect on the computed solution for Reg* > 1500. 


The upstream velocity profiles, which must be prescribed in the 
numerical computation were obtained from Cole's wall-wake law, (Ref. 38), 


u 




+ C + 




u 


+ 


where 


Wall-Wake Law 

► 

laminar sublayer 


(41) 




u r * /p 


and H(x)/< is evaluated from the condition that u = u» at y/6 = 1. In 
Eq. (41) the constants < and C are set at .41 and 5.0 respectively. This 
leaves two free parameters T w (or Cf) and 5, which must be chosen to 
completely specify the profile. Since Reg* is to be specified instead of 6 
an additional relationship is required. In Ref. 38 a relationship 
relating Reg*, 6, Cf and II(x) is given 


k ( Re^g. 65) 
8u t /j/ 


1+ II (x) 


(42) 
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Combining Eq. 41 and 42 II(x) can be eliminated and a transcendental equation 
for X = 6u t /v is obtained (with Cf and Reg* being specified constants), 

^,± lnX+ 4- (Re 65) .A (43) 

U T K X 8 K 

which can be solved for X and hence 6 by a Newton iteration procedure. In 
general, only two or three iterations are required for convergence. There- 
after <5 and Cf are substituted in Eq. (41) to obtain the required velocity 
profile. In some cases the 6 and Cf obtained from the iteration procedure 
were used as input for Musker's profile (Ref. 40). The advantage of this 
profile is that there is a smooth transition between the laminar sublayer and 
the logarithmic portion of the boundary layer. 

Streamwise locations for the upstream and downstream boundaries were 
chosen consistent with the boundary layer properties at these locations. 

This wa$ done in order to present the unsteady data as a function of ujx/u 0 , 
the reduced frequency. However, the x value so determined will not 
necessarily correspond to the x locations in Karlsson's experiment due to 
virtual origin effects which are consequences of upstream history. A 
correlation between Cf and Re x was employed, to obtain these locations. 

The correlation which assumes a l/7th power law profile for the velocity is 
given by 

C f = 0.0592 

Consequently, in order to obtain a Reg* in excess of 4000 at the downstream 
boundary, the extent of the computational domain in the streamwise direction 
was set as .594 m x _< 4.862 m (1.95 ft x 15.95 ft). 

In the normal to the wall direction, the outer boundary was chosen so 
that during the oscillatory motion the boundary layer edge would lie totally 
within the computational domain. Hence, the outer edge was placed at a 
nominal height of X 3 = 13.716 cm (.45 ft). In the normal direction either 41 
or 51 grid points were employed. For the 41 point grid a clustering 
transformation of the hyperbolic tangent type was used, with the first grid 
point being located at X 3 = .003888 cm (.000324 ft) or y + = 1.407 at the 
upstream boundary. For the 51 point calculation Oh's transformation 
(Ref. 41) was used, which relies upon an error function series for the grid 
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clustering. The first grid point was located at X 3 = .005732 cm (.000188 ft) 
or y + = 1.0377. There were 36 uniformly spaced grid points in the stream- 
wise direction yielding a mesh spacing of Axi = 12.192 cm (.40 ft) and cor- 
responds to nearly six boundary layer thicknesses as measured at the upstream 
station. 

For the quasi-steady state case in order to correspond to the 
experimental data the Reynolds number per unit length was, R e = 307,579/m 
(93,750/ft) where U» = A. 572 m/sec (15 ft/sec) and v = .0000149 sec/m 2 
(.00016 sec/ft ). For the other cases, the Reynolds number per unit length 
was Re = 358,842/m (109,375/ft) where U,*, = 5.334 m/sec (17.5 ft/sec). The 
mesh employed in the calculation consisted of 51 x 36 = 1836 grid points. 

The calculation reached a steady state in 42 time steps, defined as the point 
at which the dependent variables u and u do not change by more than e = 

10 over two successive time steps* In obtaining the steady state, the 
solution is advanced in time from some initial guess, and the time steps are 
chosen in order to hasten convergence* 

A comparison between the predicted mean velocity profile at x = 4.252 m 
(13.95 ft) corresponding to (Re<$* = 3551 and Cf = .00331) with Karlsson's 
data is shown in Fig. 7. The predicted values are in excellent agreement. 

For the unsteady calculations comparisons were made with the data at this 
same location. 

Unsteady Calculations 

The procedure for obtaining unsteady flows is similar to the steady 
solution procedure. Whereas in the steady state case, the inflow and outer 
edge boundary conditions are perscribed to be invariant in time, for the 
unsteady case the velocities at these boundaries are allowed to change in 
time. For the outer boundary this poses no difficulty since the streamwise 
velocity is known, i.e. 


u = u 0 (l + Acoscut) 


However , there is some difficulty at the upstream boundary where the stream- 
wise velocity is unknown. Several methods have been used by investigators 
(Ref. 33 - 36) to handle the specification of these upstream conditions. In 
references 34 and 35, the authors begin their calculation at x\ = 0 as a 
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laminar flow and arbitrarily specify a location for instantaneous transition, 
to fully turbulent flow. Although such a procedure can be employed when the 
governing equaions are cast in boundary layer coordinates (e.g., n=y/ /2x) in 
order to remove the singularity at x\ = 0, it necessitates the consideration 
of transition which we have attempted to avoid. Others, (Ref. 33 and 36) 
have specified quasi-steady upstream conditions, in which Cf and 6 are 
varied (the parameters in the Cole's velocity profile) to correspond to a 
steady state that would exist with the same edge velocity at that instant in 
time. This method also has drawbacks, since the upstream profile is void of 
any phase shift effects, and this is incompatible with the rest of the flow 
downstream of that location. 

In the present procedure the upstream profile is still given by the 
Cole's wall-wake velocity profile, but now Cf and 6* are allowed to be 
periodic functions of time; l.e.. 



8* = 8*( I + AS* cos (tut + 7 r + ) 

where Cf Q and 6* 0 are the mean (time averaged) skin function and dis- 
placement thickness, Cf^ and 6*^ are their respective amplitudes of 
oscillation and <j> C f and <f>6* their respect phase shifts. This procedure 
introduces additional unknowns. For the mean quantities Cf 0 and 6* 0 , the 
steady state values are used, while the other four quantities are determined 
by analyzing the solution behavior near the upstream boundary. For instance 
the variation of 4> c ^ as a function of the reduced frequency 0ix/u o is 
shown in Fig. 8. As can be seen, the quasi-steady upstream profile is not 
consistent with the interior, while the present procedure gives more 
realistic behavior leading to a smoother transition from the upstream to the 
interior. Although wx/u 0 may not be the unique scaling parameter for 
turbulent flow, it is nevertheless used in the present approximate analysis 
which is now briefly described. First, the parameters in Eq. (44) were 
assumed for the low frequency case, to = .33 Hz, and a solution was computed. 
Thereafter, the distributions of cf^, $cf, 6*i> an< * $6* were obtained 
as a function of u)x/u Q in the initial calculation. For another frequency, 

^2 for example, the corresponding values of Cf 1# 4> C f> 6* j and 4>$* 

can be determined at the new value of 0 ) 2 X*/u o , from a curve such as shown 

in Fig. 8, where x* is the location of the upstream boundary. 
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For all the reported computations a time step corresponding to 10° of 
the period was chosen. The calculation was run through 2 periods (72 time 
steps) and the results were obtained for the second period. Several 
calculations were run through 3 periods and the solutions did not vary from 
the second period results. Also, several calculations were run in 5° time 
step increments and again the results were almost indistinguishable from 
those obtained with 10° time steps, thus indicating accurate temporal 
resolution. The run time for a 36 x 51 grid is 44.91 sec on the CRAY-1 for 
three cycles (108 time steps) and corresponds to .00024 sec/grid point/time 
step. 


In presenting the results, the skin friction, Cf, the displacement 
thickness, 6* and the streamwise velocity profile at the measuring station 
were Fourier decomposed into their harmonic components. 

For a function f(*r), the Fourier Series becomes: 


0o oo * 

f(r) = — +1 |a n coscunt+ b n sinajnt} 
Z n=l 


where 


t, + T 

°o»TF /«€><!£ 


U) 

°n = tt 


t . +T 


J 1 (£)cosGjn£d£ 


(45) 


t,+T 


b n ™ f f(Osinu>n£d£ 


and where ti is the time at the start of the integration, and the period T = 
2tt/(d. 

Although the code allows for the determination of any number of 
harmonics, only the first two were obtained. In Karlsson's experiment, the 
mean (time averaged) velocity and the components of the first harmonic 
were measured. Therefore, comparing to the Fourier series representation the 
mean velocity is a Q /2, the in phase component of the first harmonic is a^ 
and the out of phase component of the first harmonic is -b \ (to correspond to 
Karlsson's notation). For the evaluation of the Fourier coefficients 
Simpson's integration scheme was used and all data points for the second 
cycle were sampled. 

Results of the computations and comparisons with the data of Karlsson 
are shown in Figs. 9-22. In the first set, Figs. 9-13, the time averaged 
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mean velocity profiles of the location corresponding to Karlsson's streamwise 
measuring station are shown for frequencies w = .33, 1.0, 1.33, 2.0, 4.0 and 
7.7 Hz with amplitudes of oscillation A = .2020, .195, .1763, .1358 and .1273 
respectively. As can be seen in these figures, there is excellent agreement 
between the computed profiles and the data. Since there is little difference 
between the mean velocity profiles and the quasi-steady profiles, unsteady 
effects cannot be inferred from an inspection of these figures. Therefore, 
the in phase and out of phase components of the first harmonic of the 
oscillations must be considered; corresponding to normalized values of ax and 
-bx in Eq. 45. These are given in the next set. Figs. 14-18. 

As can be seen in these figures, there are considerable unsteady 
effects present with the predominant effects appearing at the lower 
frequencies. At the highest frequency calculated, 7.7 Hz, the unsteady 
effects are relegated to a thin region near the wall reminiscent of the shear 
wave solutions of Lighthill (Ref. 37) for laminar flows. Note that the 
boundary layer does not oscillate as a unit at the excited frequency but 
rather different portions oscillate at varied frequencies. This is further 
exemplified by the characteristic overshoots and undershoots in the in phase 
and out of phase components respectively that were obtained for all cases 
run. For m = .33 Hz there is excellent quantitative agreement with data for 
both components, while for the other frequencies considered the correct 
qualitative behavior was obtained but the exact magnitudes and locations of 
the peaks were not recovered. It should be noted that the present results 
show quite good agreement with Karlsson's data even though a simple 
turbulence model was. used. In particular, the characteristic in phase and 
out of phase velocity profiles were qualitatively predicted and 
quantitatively these profiles were in agreement with the data. The existing 
descrepancies , however, could be due to turbulence modeling, upstream 
conditions or experimental error. It should also be noted that the present 
in phase and out of phase prediction show better agreement with Karlsson's 
data than was obtained by other researchers using a variety of turbulence 
models (e.g. Refs. 42 and 43). 

The next two figures 19 and 20 show the amplitude and phase angle 
variations of the displacement thickness plotted against the reduced 
frequency u)x/u Q . Also shown in these figures are the data of Karlsson 
(interpreted by Telionis with his error estimates indicated by vertical bars) 
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as well as the present predictions for oi « *33, 1.0, 2.0, 4.0 and 7.7 Hz. 

Good agreement is obtained between the present results and the data as far as 
amplitude is concerned. In addition, these results compare favorably with 
those obtained by McCroskey (Ref. 33). 

However, the phase angle of the skin friction (phase lead) as a 
function of reduced frequency (cf. Fig. 21) shows an apparent lack of 
agreement between the present predictions and the data. Such a discrepancy 
was also noted in the results of other analyses as shown previously in 
Fig. 4. A better understanding of these results can be obtained by viewing 
Fig. 21 which presents the calcuated c)> c f versus <j4x/u q for each of the 
five runs made under the present effort. Each curve represents the results 
for a different oscillation frequency, 

As can be seen in Fig. 21 for each of the frequencies the phase angle 
approaches a different constant asymptotic value that is increasing with 
frequency. Each case taken alone resembles the behavior obtained for laminar 
flow. Whereas in laminar flow the unique similarity variable is the reduced 
frequency, a)x/u 0 , for turbulent flow it may not be. Furthermore, the phase 
angle is likely to be a function of the Reynolds number as well. Therefore, 
when viewed in this manner, the present predictions would have the 
appropriate physical behavior since they mimic the results obtained for 
laminar flow. Nevertheless, the lack of correlation with the data and the 
correspondence of the various calculations still needs to be addressed. 

The experimental data as well as the other computations in Fig. 4 tend 
to lie on curves that are monitonical ly increasing with reduced frequency. 
Recall that the experimental data correspond to different frequencies at one 
streamwise location, and not to the phase angle distribution along the plate 
at a given frequency. If one would connect the points associated with the 
measuring station for each of the present calculations one would obtain a 
curve resembling for instance that of McCroskey (Ref. 33) (cf. Fig. 4). 

The curves obtained by other investigators appear to be based on data points 
for different values of to and may not indicate the appropriate physical 
behavior, nor for that matter they' may not contradict the present results. 
Finally, with regard to the displacement of the various curves, the 
discrepancy may reside in the virtual origin associated with turbulent 
flows. Since the streamwise distance is not a characteristic length, the 
distances in various predictions would not correspond precisely with those of 
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Karlsson. By taking this into account, and shifting the present data points 
so that for the w = . 33 Hz case the reduced frequency of the computation is 
set equal to that of the experimental data the new curve (stars) in Fig. 22 
is obtained. 


Three-Dimensional Calculation 

The three-dimensional unsteady flow requires additional attention 
compared to its two-dimensional counterpart, ir particular with respect to 
boundary condition specification. For two-dimensional flows at the upstream 
inflow boundary (which is a line) the streamwise velocity profile must be 
specified as a function of time. For the three-dimensional flow, there are 
two inflow boundaries (planes) (cf. Fig. 2) where streamwise as well as 
spanwise velocity profiles must be specified. Since the intent of this 
calculation is to demonstrate the capability of computational procedure for 
three-dimensional unsteady flows, the boundary conditions at these boundaries 
were obtained from the two-dimensional results. The description of the 
actual steps in the calculation are now given. 

First, the analogous two-dimensional steady and unsteady calculations 
were completed. The normal and streamwise velocity profiles were then stored 
for each time step during the second cycle of the oscillating flow. This 
corresponds to the flow along the diagonal of the square computational 
domain. For the two-dimensional problem the extent of the computational 
domain was .960 m xi 4.526 m (3.15 ft < xx < 14.85 ft) and 
0 < *3 < .122 m (0 C x 3 £ .4 ft) and a 27 by 41 grid was employed. In the 
streamwise direction, a uniform grid was employed with Axx = .137 m (.45 ft) 
while normal to the wall a hyperbolic transformation was used to cluster 
points near the wall. The setting of the upstream profile was discussed 
previously. In anticipation of the three-dimensional problem to be done the 
upstream xx location was located at xx = ,960 m (3.15 ft) rather than .594 m 
(1.95 ft). Therefore, the Cole's velocity profile required different values 
of Reg* and Cf which were set to 1200 and .0043, respectively. Steady 
state was reached in 38 time steps with a total running time of 9.545 seconds 
or .00023 sec/grid point/time step. Thereafter the unsteady case was run for 
w = . 33Hz and with an amplitude of oscillation, A = .202. The modified 
upstream condition procedure was used, with the following parameters set 
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6*x = .28 Cfl = 1.82 

<f>6* = 24“ 4>cf = 4° 

The calculation was run for two periods with a time step corresponding to 10° 
and the Fourier decomposition was computed on the second cycle. Results of 
the calculation are shown in Figs. 9 and 14. The mean flow and the in-phase 
and out-of-phase components are in excellent agreement with the data. The 
phase shifts angles and amplitudes show qualitatively the same behavior as 
the calculations of McCroskey and Phillipe (Ref. 33). 

Having obtained the two-dimensional unsteady results, a three- 
dimensional flow field calculation can be made. The first task is to 
determine a starting flow field for the three-dimensional calculation. Since 
a periodic solution is sought, the flow field at the beginning of a cycle 

will not necessarily correspond to the steady state flow field, but rather to 

the cyclical flow. In order to reduce total computational time, the initial 
flow field used to run the 3-D calculation was that which was obtained at the 
beginning of the second cycle of the two-dimensional calculation. Recall 
that the two-dimensional calculation corresponds to the flow along the 
diagonal of the computational domain in the 5 coordinate direction. Hence, 
in order to obtain the initial flow field the two-dimensional data was 
reflected directly onto the grid by decomposing the 5 direction velocity into 
its xi and X 2 components. This was done in order not to interpolate the 
two-dimensional calculation onto the three-dimensional grid. This resulted 
in a 14 x 14 grid in the xi~X 2 plane and a grid spacing of Axi=Ax2= .194 m 
(.636 ft). Since this spacing is larger than that used for the 
two-dimensional calculation, there is a question to what effect that 
discretization error would have on the solution. Therefore, another coarser 
two-dimensional calculation was performed on a 20 x 41 grid which 
corresponded to the same Axi as the three-dimensional calculation, and the 
results obtained were very close to the fine mesh calculation. 

Once the initial flow field was obtained, the three-dimensional 
calculation was run, also with a 10° time step, but for only one cycle. 

Since the initial state is obtained from a periodic solution there is no need 
to march through one preliminary cycle in order to eliminate transient 
non-periodic effects. The inflow boundary conditions as stated previously 
were the two-dimensional unsteady results, decomposed into their xi and X 2 
components, and applied along boundaries 1 and 2 in Fig. 2. 
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The results for the mean flow and the in phase and out of phase 
components of the first harmonic of the velocity were indist ingushable from 
the two-dimensional results (cf. Figs. 9 and 14). This indicates that the 
three-dimensional procedure can give accurate results comparable to the 
two-dimensional calculation. However, further validation was obtained by 
comparing the three-dimensional derived boundary layer properties with the 
two-dimensional results along a ? = constant line (perpendicular to the 
diagonal). Along such a line the values should be equal to the two- 
dimensional results. In Table 1 these quantities are presented as a function 
of plane location corresponding to points along the line perpendicular to the 
diagonal in the adjacent sketch. 

In the first column, the two integers signify the X2 and xi grid 
locations in that order. Also given in the table are the 27 x 41 fine and 
20 x 41 coarse, two-dimensional calculations. As can be seen, the values 
along the diagonal are relatively constant and the maximum relative error 
betwen the three-dimensional and the fine two-dimensional results are less 
than 0.6%. The minimal variation of these derived values along the 
£ = constant line indicates that the three-dimensional calculation procedure 
recovered the two-dimensional results, hence, further validating the 
technique. Since the solution is ostensibly two-dimensional in the 
appropriate frame of reference, the physical behavior corresponds to that 
which was discussed previously and thus the reader is referred to that 
section. 

As a final note, the reader is referred to Table 2 where run times 
(total and per grid point per time step) are given for both two-dimensional 
and three-dimensional results. These run times are for the unvectorized 
code. If past experience with vector izat ing other fluid mechanics codes 
holds for this code then one can expect a reduction of up to a factor of 
five in run time. Hence, it is conceivable that three-dimensional viscous 
calculations of the type considered in this report can be obtained in less 
than a minute of CPU time on the CRAY-1 computer. 
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CONCLUSIONS 


In the Phase I effort, all major objecties have been satisfied. The 
three-dimensional solution procedure for the approximate form of the 
Navier-Stokes equation was exercised in the two- and three-dimensional modes 
to compute the unsteady turbulent boundary layer on a flat plate 
corresponding to the data of Karlsson. New upstream boundary conditions were 
developed that yielded more realistic solutions for the interior in the 
vicinity of the upstream boundary. Comparisons of the computaion employing 
these boundary conditions with the data indicate that both qualitative and 
quantitative agreement was obtained for the mean velocity and the in phase 
and out of phase components of the first harmonic of the velocity. 

In addition, the calculation gave results for the skin friction phase 
angle that had plausible physical behavior for large distances downstream of 
the inflow boundary. Finally, an explanation was suggested to resolve the 
discrepancy with regard to previously reported calculations of the skin 
friction phase angle. 

For the three-dimensional case, the two-dimensional data of Karlsson 
was considered, but skewed at 45° by a coordinate rotation. The results of 
the calculations were in excellent agreement with the data and the 
two-dimensional computations, thereby validating the three-dimensional 
procedure . 
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Figure 1 - Effect of Coordinate System on Grid Distribution 



Figure 2 - Three-Dimensional Computational Domain 
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Figure 4 - Phase Lead Angle of Skin Friction as a Function of Reduced Frequency 











Figure 6 - Skin Friction Coefficient as a Function of 
Displacement Thickness Reynolds Number 
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Figure 7 - Quasi-Steady Velocity Profile 
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Figure 8 - Effect of Upstream Boundary Condition 















Figure 14 - In Phase and Out of Phase Velocity Components 














Figure 18 - In Phase and Out of Phase Velocity Components 
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Figure 19 - Amplitude of Displacement Thickness Fluctuations 
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Figure 20 - Phase Angle of Displacement Thickness Fluctuations 
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Figure 21 - Phase Angle of Skin Friction Fluctuations 





Figure 22 - Phase Angle of Skin Friction Fluctuations - Effect of Virtual Origin 
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Table 1 - Comparison of Three-Dimensional Results 


















ON 



GRID 

TIME STEPS 

TOTAL 

PER GRID POINT /TIME STEP 

2-D 

36 X51 
(1836) 

108 

44.91 

.00024 

3-D 

14 X 14 X 41 
(8036) 

36 

185.35 

.00064 















1. Report No. 

NASA CR-172368 


2. Government Accession No. 


3. Recipient's Catalog No. 


4 . Title and Subtitle 

THREE-DIMENSIONAL UNSTEADY VISCOUS FLOW ANALYSIS OVER 
AIRFOIL SECTIONS 


5. Report Date 

June 1984 


6. Performing Organization Code 


7. Author(i) 

B. C. Weinberg and S. J. Shamroth 


9. Performing Organization Name and Address 

Scientific Research Associates, Inc. 
P.0. Box 498 
Glastonbury, CT 06033 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 


8. Performing Organization Report No. 

R84-900027F 


10. Work Unit No. 


11. Contract or Grant No. 

NAS1-17573 


13. Type of Report and Period Covered 
Contractor report 


14. Sponsoring Agency Code 


15. Supplementary Notes 


Langley Technical Representative: John W. Edwards 

Final Report 

16. Abstract 

A three-dimensional solution procedure for the approximate form of the Navier 
Stokes equation was exercised in the two- and three-dimensional modes to compute 
the unsteady turbulent boundary layer on a flat plate corresponding to the data of 
Karlsson. The procedure is based on the use of a consistently split Linearized 
Block Implicit technique in conjunction with a QR operator scheme. New time- 
dependent upstream boundary conditions were developed that yielded realistic solu- 
tions for the interior in the vicinity of the upstream boundary. Comparisons of 
the computation employing these boundary conditions with the data indicate that 
both qualitative and quantitative agreement was obtained for the mean velocity and 
the in phase and out of phase components of the first harmonic of the velocity. 

In addition, the calculation gave results for the skin friction phase angle that 
had expected physical behavior for large distances downstream of the inflow 
boundary. For the three-dimensional case, the two-dimensional data of Karlsson 
was considered, but in a coordinate system skewed at 45° to the free stream 
direction. The results of the calculations were in excellent agreement with the 
data and the two-dimensional computations. 


17. Key Words (Suggested by Author(s)) 18. Distribution Statement 

Three-Dimensional 

Unsteady Unclassified - Unlimited 

Turbulent 



«-}« 


For sale by the National Technical Information Service. Springfield. Virginia 22161 



















